clear; close all; clc
warning off;
 addpath('../LibForSerial/LibrarySLP','../LibForSerial/LibraryLeSage/','../LibForSerial/LibSoderlind/','../LibForSerial/LibVAR/');

%----------- Load shocks using two identification schemes ----------------%

NN=9;
Tmax = 117; %86, 110 for France (6), 116 for Japan (2), 117 for CHE (10)
 ydk = NaN(Tmax,NN);
 zdk = NaN(Tmax,NN,10); % w/  employment
nn=0;

% Settings
%WHICH MACRO UNCERTAINTY INDEX? Baker, Bloom & Davis (2016) (BBD) OR JLN
muindex=0; % 1 =Macro U ; 0 = BBD
%COUNTRY LABELS:
clabel ={'USA','Japan','Germany','Italy','UK','France','Canada','Spain',...
    'Sweden','Switzerland','Netherlands'};

for country=[5 7 6 3 2 4 8 9 10 11]
    strcat('Country = ', clabel{country})
%-------------------------------------------------------------------------%
LogVsLevel = 1; 
  NumberConsecutiveshocks = 2; % 1 previous shocks that is  positive, and the 2nd one is positive too 

cd ./DataShocks/
LoadShocksMandQpositiveEPUrevisedMAC; 
cd ../ 
%-------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data = dataMatrix;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA TRANSFORMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % data(:,[3:10]) = log(data(:,[3:10]))*100;
%--------------------------------------%
clear startsample endsample 

T = size(data,2);
%P  =  1;
%P  =  4; % including in the LP regression P lags of all variables in the system (3 lags in original BBD QJE; here 6 for consistency througout)
 P  =  2; % 
                            
 plottype = 'Output';      ii=3; % 
 
 %------------ IRF of GDP growth to a shock variable ----------------------%
 y  = data(:,ii);    %endogenous variable

 %-------------------------------------------------------------------------%
  Remove_Large_Shocks = 0; 
 % G7 + Spain + Sweden + Switzerland + Netherlands
 %  1. USA, 
 %  2. JP , Insignificant, but OK
 %  3. GER, 
 %  4. IT, 
 %  5. UK, 
 %  6. FRA, 
 %  7. CN,  
 %  8. ESP, 
 %  9. SWE, TO BE FIXED
 % 10. CHE, Insignificant, but OK
 % 11. NL , 
%  if country == 5 || country == 7
%      Remove_Large_Shocks = 1;
%  end
    x  = uEPU;  
    IndSeqShocks = IndPreviousNShocksEPU;
    IndLargeShocks = uEPU>2;
    
    LargeShocksToBeRemoved = IndSeqShocks & IndLargeShocks;  

 
  %fprintf('Number of %d consecutive shocks to be removed  %d\n', NumberConsecutiveshocks, sum(LargeShocksToBeRemoved));
   % [x IndSeqShocks IndLargeShocks]
   if Remove_Large_Shocks; IndSeqShocks = IndSeqShocks - LargeShocksToBeRemoved; end
  %fprintf('Number of events with %d consecutive shocks is %d \n',NumberConsecutiveshocks,sum(IndSeqShocks));
  
  %/////////////////////////////////////////////////////////////////////////%
  %-------------------------------------------------------------------------%
   xs = x.*IndSeqShocks; %shock variable times states
  %-------------------------------------------------------------------------%
 % w  = [lagmatrix( data(:,unique([ii 2       ])  ), 1:P )                       ]; %control variables (contemporaneous vars, lagged vars)
 % w  = [lagmatrix( data(:,unique([ii 2       ])  ), 1:P ) lagmatrix(logEPU_rescaled,1:P )]; % w/o employment
   w  = [lagmatrix( data(:,unique([ii 2   4   ])  ), 1:P ) lagmatrix(logEPU_rescaled,1:P )]; % Same controls as in the 1st step VAR

OKforReg = isfinite(x) & isfinite(y) & isfinite(xs);
for cc=1:size(w,2)
    OKforReg = OKforReg & isfinite(w(:,cc)) ;
end  
  x = x(OKforReg,:); 
  y = y(OKforReg,:)*100; 

  w = w(OKforReg,:);
  xs=xs(OKforReg,:);
  
  
% Do a simple Panel AR ... 
nn=nn+1;
regressors = [x xs w];%[x w];
if country ~= 5 
    ydk(end-size(y,1)+1:end,nn)   = nanstd(ydk(:,1))*(y-nanmean(y))/std(y) + nanmean(ydk(:,1));
    regressors(:,1)               = nanstd(squeeze(zdk(:,1,1)))*x /std(x);
    regressors(:,2)               = nanstd(squeeze(zdk(:,1,2)))*xs/std(xs);
    for rrr=1:8%8 w/ employment, 6 w/o employment
    regressors(:,rrr+2) = nanstd(squeeze(zdk(:,1,rrr+2)))   *(w(:,rrr)-nanmean(w(:,rrr)))/std(w(:,rrr)) + nanmean(zdk(:,1,rrr+2));
    end
else
    ydk(end-size(y,1)+1:end,nn)   = y;
end
zdk(end-size(y,1)+1:end,nn,:) = regressors;

end

%% IRF
h1 = 0 ;  %start LP at h=0 or h=1 (h=1 if impose no contemporaneous impact)
H  = 4*3; %#horizon of IR

%   BW = norminv(0.975); % 95% confidence level --> 1.96
    BW = norminv(0.95) ; % 90% confidence level --> 1.64
    
maxhorz = H; 
NWlags  = 12;
regLin       = NaN(maxhorz,1);
regLinState  = NaN(maxhorz,1);
confidenceLin      = NaN(maxhorz,2);
confidenceLinState = NaN(maxhorz,2);
[A,stdDK,stdW,CovDK,yhat,R2,CovW,CovDKj,B] = HszDk5cPs(ydk,ones(length(ydk),1),zdk,1,NWlags);
regLin(1)      = A(1);
regLinState(1) = A(2);
confidenceLinState(1,:)=[A(2)-(B(2)*BW),A(2)+(B(2)*BW)];
for hh=2:maxhorz+1
    ydk_h = ydk(hh:end,:);
    zdk_  = zdk(1:end-hh+1,:,:);
    ydk_h(:,sum(~isfinite(ydk_h))==size(ydk_h,1))=[];
    zdk_( :,sum(~isfinite(ydk_h))==size(ydk_h,1))=[];
    [A,stdDK,stdW,CovDK,yhat,R2,CovW,CovDKj,B] = HszDk5cPs(ydk_h,ones(length(ydk_h),1),zdk_,1,NWlags);
   
           regLin(hh,1)= A(1);
    confidenceLin(hh,:)     =[A(1)-(B(1)*BW),A(1)+(B(1)*BW)];
           regLinState(hh,1)= A(2);
    confidenceLinState(hh,:)=[A(2)-(B(2)*BW),A(2)+(B(2)*BW)];
end



ScaleFactor = 0.5; % To normalize maximum drop to 1% as in Redl
%-------------------------------------------------------------------------%
%-------------------------------------------------------------------------%   
% For compatibility with previous code
%-------------------------------------------------------------------------%
slp_linear = regLin;
slp_sd.IR  = regLin;
slp_sd.IRstate=regLinState;
slp_sd.Interval_down_state = confidenceLinState(:,1);
slp_sd.Interval_up_state   = confidenceLinState(:,2);

subplot(1,2,1)
  plot( 0:H ,ScaleFactor*(slp_linear), 'k', 0:H, ScaleFactor*(slp_sd.IR+0.25*slp_sd.IRstate), 'b--', 0:H, ScaleFactor*(slp_sd.IR),'k--', 'LineWidth', 1.5 ); 
  title(['IR_{slp} when there are ' num2str(NumberConsecutiveshocks) ' Consecutive Positive Shocks']); %legend(         's_t=+1 (Previous shock(s)=1)','s_t=0','Location','Best')
hold on; plot(0:H , zeros(H+1,1) , '-k' , 'LineWidth' , 2 ); grid; xlim([0 H]);                         legend('Linear','s_t=+1 (Previous shock(s)=1)','s_t=0','Location','Best')
if ii==2 
    ylim([-2.8 1.6]);set(gca,'YTick',[-2.8:0.5:0.6],'FontSize',12);ylabel('Percent','FontSize',12);
% elseif ii == 4
%     ylim([-1.5 2]);  set(gca,'YTick',[-1.5:0.5:2]  ,'FontSize',12);ylabel('Percent','FontSize',12)
end
xlim([0      H]);set(gca,'XTick',[0 3:3:H],'FontSize',12);

subplot(1,2,2)
grpyat=[(0:1:H)', ScaleFactor*0.25*slp_sd.Interval_up_state; (H:-1:0)' ScaleFactor*0.25*slp_sd.Interval_down_state(end:-1:1)]; patch(grpyat(:,1), grpyat(:,2), [0.7 0.7 0.7],'edgecolor', [0.7 0.7 0.7]); hold on
plot(0  :H , zeros(H+1,1), 'k-'); hold on 
plot(    0:1:H,   ScaleFactor*0.25*slp_sd.IRstate, 'b--','LineWidth', 1.5); hold on
title('State multiplier of State-dependent IR to uncertainty shock')
grid;
xlim([0      H]);set(gca,'XTick',[0 3:3:H ],'FontSize',12);